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Abstract 

The thin plate spHne is a popular tool for the interpolation and smoothing of 
scattered data. In this paper we propose a novel stabilized mixed finite element 
method for the discretization of thin plate splines. The mixed formulation is obtained 
by introducing the gradient of the smoother as an additional unknown. Working with 
a pair of bases for the gradient of the smoother and the Lagrange multiplier which 
forms a biorthogonal system, we can easily eliminate these two variables (gradient of 
the smoother and Lagrange multiplier) leading to a positive definite formulation. The 
optimal a priori estimate is proved by using a superconvergence property of a gradient 
recovery operator. 
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1 Introduction 

We propose a new finite element approach for the discretization of the the thin plate spline 
[18,33], which is one of the most popular approach in scattered data fitting. Scattered 
data fitting problems occur in many applications such as data mining, reconstruction of 
geometric models, image processing, parameter estimation, optic flow, etc., see [6,21,34]. 

Let Q, C M'* with d G {2, 3} be a closed and bounded region with polygonal or polyhe- 
dral boundary. In the following, we use the standard notation for the norm and semi-norm 
of Sobolev spaces [12]. Given a set Q = {xjj^j^ of scattered points in Q, and a function 
r on ^ with Zi = r(xj) for i = I,-- - ,N, the thin plate spline is a smooth function 
u-.n^m [18,33] such that 

min F(u), (1) 
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where ^ 

F{u) = f^i^i^i) - ^i? + « / E f^l (^'^)' (2) 
1^1 1^2 ^""^ 

u = {vi, ■ ■ ■ , Vd) G Ng is a multi-index, \v\ = Yli=i a is a positive constant. 

A conventional approach is to use radial basis functions to approximate the space 
H'^{Q) in ([1]), which leads to a dense system matrix. The solution of such a system is very 
expensive when a large data set has to be modelled. In this paper we propose an efficient 
discretization technique for the minimization of the functional ([T]). The basic idea of a 
finite element method is to minimize the functional F given by ([2]) over a finite-dimensional 
function space. If we want to discretize the minimization problem using a conforming 
approach, we need to construct a discrete finite element space which is a subset of the 
Sobolev space H'^{Q). Construction of such a finite element space is expensive [12, 16]. 
The class of standard non-conforming finite elements [10, 16] provides a more efficient 
discretization than the conforming approach. However, their implementation requires a 
more complicated data structure, and a suitably constructed mixed formulation provides 
a more efficient and fiexible discretization than the non-conforming approach. Therefore, 
following a similar approach as in [2, 15, 22, 29], we modify the original minimization 
problem ([1]) so that the minimization is done over the Sobolev space H^{Q) rather than 
over the Sobolev space H'^{Q), and the formulation allows an efficient mixed finite element 
discretization. A similar idea has been exploited in [16, 17, 19, 28] for the solution of 
biharmonic equation with simply supported and clamped boundary condition. 

The rest of the paper is organized as follows. In the remainder of this section, we fix 
some notations and introduce an alternative formulation. The next section is devoted to 
the finite element formulation of the problem. We put the problem in a framework of 
an appropriate saddle point problem, and discuss its algebraic structure. The algebraic 
structure of the problem motivates us to use a pair of finite element bases for the gradient 
of the smoother and the Lagrange multiplier which forms a biorthogonal system. Section 
[3] is devoted to the analysis of the discrete problem. Eliminating the gradient and the 
Lagrange multiplier, we get a positive definite formulation of the saddle point problem for 
which we prove the existence of a unique solution. The final part of this section shows the 
optimal convergence of our finite element solution to the continuous solution. We conclude 
the paper with a summary. 

Let the Sobolev space H^{Q) x [i7^(r2)]'^ be denoted by V, and for two matrix- valued 
functions a : J7 ^ M^^"' and /3 : 1^ ^ M'^^"', the Sobolev inner product be defined as 

d d 

: = E E j 1 ) gfc (H) ; 

i=i j=i 

where {oi)ij = aij, {l3)ij = with aij,Pij G H^{^1), and the norm || • is induced 

from this inner product. For A: = 0, an equivalent notation 



(a,/3)L2(n) := 




2 



for the L^-inner product will be used and the L^-norm || • ||i2(f^') is induced by this inner 
product. 

A new formulation of the functional F in ([1]) is obtained by introducing an auxiliary 
variable cr = Vu such that the minimization problem ([1]) is rewritten as [15, 22] 

min G(u, cr) , (3) 
{u,a-)&v 
cr=Vu 

where 

N 
i=l 



2 Finite element approximation 

Let Tfi he a quasi-uniform partition of the domain 0, in d-simplices or d-parallelotopes 
having the mesh-size h. Let T be a reference d-simplex or d-cube, where a reference 
d-simplex is defined as 

d 

f := {{xi, ■ ■■ ,Xd) eW^ : xi > 0, • • • ,Xd> and < 1}, 

i=l 

and a reference d-cube T := (—1,1)'^. 

The finite element space is defined by affine maps Ft from a reference d-cube or d- 
simplex T to a d-parallelotope or a d-simplex T G T^. Let Q{T) be the space of multilinear 
polynomials in T if T is a reference d-cube or the space of linear polynomials in T if T is a 
reference d-simplex. Then the finite element space based on the mesh is defined as the 
space of continuous functions whose restrictions to an element T are obtained by maps of 
multilinear or linear functions from the reference element; that is, 

Sh := {vh G H^n) : Vhk = Vh o Vh G Q(f ), T e , (4) 

see [10,12,16]. 

Let Mh C L^(il) be another finite element space satisfying the following assumptions. 
Assumption 1. U\(i) dimAf^ = dim 5/1. 

U^ii) There is a constant /? > independent of the triangulation such that 

Uh\\L^m<P sup kl^^h^^ ^h^Sh. (5) 



U^iii) The space has the approximation property: 



inf U - AhllL^CQ) ^ Ch\<p\H^n), ^ e H\n). (6) 
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To obtain the discrete form of the minimization problem we introduce a finite 
element space V^, which is a discrete counterpart of V as Vh = Sh x [Sh]'^- Replacing the 
space y in ([3]) by our discrete space Vh, our discrete problem is to find 



N 

min 

{uh,CT h)&Vh '.^-^ 



Y^{uh{xi) - Zif + a\\V(Th\\l2(^^^ (7) 
subject to 

{o-h,Th)L^(n) = {^Uh,Th)L\n), Th G [Mhf. (8) 
If we modify the constraint ([8]) to 

we obtain the finite element thin plate spline presented in [2,30]. There are two drawbacks 
of the finite element thin plate spline presented in [2,30]. The first one being the saddle 
point structure of the system matrix arising from the discretization which is difficult to 
solve. The second drawback is that it does not necessarily converge to the standard 
thin plate spline although it has similar smoothing properties as the standard thin plate 
spline [30]. 

Now we introduce an equivalent saddle point formulation of our approach, which can be 
shown to be equivalent to the minimization problem ([7]) by using the ideas in [14,16]. We 
denote the vector of function values of u G C^{^) at the measurement points xi , X2 ,■• • ,X7v 
by G M^, i.e., 

Pu = (n(xi),n(x2), • • • ,n(xAr))'^. 

Introducing a Lagrange multiplier unknown cf)^^, the variational saddle point formulation 
of the minimization problem ([7]) is to find {{uh, (Th),4>h) ^ x [MhY so that 



(9) 



Mi'^h,o-h),{vh,rh)) + B{(t)f^,{vh,Th)) = fivh), ivh,Th) G Vh, 
B{il,h,{uh,<Th)) =0, G [Mh]\ 

where bilinear forms j4(-,-), B{-,-) and /(■) are given by 

A{{uh,(Th),{vh,Th)) = {Puh)'^Pvh + a / V(Th -.VThdx, 

Bi'iph7ivh,rh)) = / Th-il)f^dx- / Vu/, • i/>;j(ix, and /(u/,) = (Pf/,)'^z. 
JQ Jn 

We recall that the mixed formulation of our problem is closely related to the mixed formu- 
lation of Mindlin-Reissner plate [4,5,9,14], and hence we use some of the ideas presented 
in [4, 14] to analyze our problem. The existence and uniqueness of the solution of the 
saddle point problem Q is performed by using the theory presented in [4, 14]. The main 
difficulty here as well as in the context of Mindlin-Reissner plate is that the bilinear form 
A{-,-) is not elliptic on the whole space Vh- However, it would be sufficient that the bilinear 
from A{-,-) is elliptic on the space KerBh defined as 



KerBh := |(vh,Th) G Vh : ^(rt - Vvh) • ^hdx = 0, G [Mh]'^| 



(10) 
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For Sh as defined by ^ and Mh satisfying Assumptions [II|i)-miii), we cannot obtain 
coercivity of A{-,-) even on the space KerBh- This gives us a motivation to modify the 
bilinear form A(-,-) consistently by adding a stabilization term so that we obtain the 
ellipticity on the space KerBh. The modification of the bilinear form A{-,-) is done as 
suggested by Arnold and Brezzi [4] for the Mindlin-Reissner plate so that our discrete 
saddle point problem is to find {{uh, crh), cj)f^) G V/i x [Mh]'^ such that 



A{{uh,CTh), {Vh,Th)) 



where the bilinear form A{-,-) is defined as 



f{Vh), 

0, 



{Vh,Th) 



G 
G 



(11) 



A{{uh,(Th), {vh,Th)) 



{PuhfPvh + a V(Th ■■ VTh dy. + r / {an - Vuh) ■ {tk - Vvh) dx 
Jn Jn 



with r > being a parameter. Since the stabilization term is consistent, the parameter r > 
can be arbitrary in principle. By choosing an appropriate parameter, the stabilization 
can, in addition, accelerate the solver as in an augmented Lagrangian formulation [8]. 
Since we do not focus on this aspect of the problem, we simply put r = 1 in the rest of 
the paper. After putting r = 1, we have 



A{{uh,(Th), {Vh,Th)) 



A{{uh,(Th),{vh,Th)) + / {(Th -Vuh) ■ {Th -Vvh)dx. 

Jn 



Here our interest is to eliminate the degree of freedom corresponding to cTh and and 
arrive at a formulation only depending on Uh- This will dramatically reduce the size of the 
system matrix, and the system matrix after elimination of these variables will be positive 
definite. For the solution of the reduced system, one can thus use very efficient numerical 
techniques. Therefore, we closely look at the algebraic formulation of the problem. In 
the following, we use the same notation for the vector representation of the solution and 
the solutions as elements in Sh, [Sh]'^ and [M/i]*^. Let R, A, B, W, K, D and M be the 
matrices associated with the bilinear forms {Puh)^ Pv^, Jq ^o^h '■ Vt/j dx, Vu^ ■ ijjh fix, 

Vuh • Tfi dx, ■ Vv/i dx, J"^ (Th ■ ipfi dx and o-^ • dx, respectively. The matrix 

D associated with the bilinear form ■ il'h often called a Gram matrix. In case 

of the saddle point formulation, u^, cth and are three independent unknowns. Letting 
the test functions Th and to be zero subsequently in the first equation of (jlip . we have 



{PuhfPvh - Vvh ■4>hd^- In^^h - ^Uh) ■ Vvh dx 

" !n "^^h ■ dx + !^4>h- '^h fix + f^{(Th - VUh) ■ Th dx 



f{Vh), 

0, 



Th G L^h 



Oh 



Then the algebraic formulation of the saddle point problem pi|) can be written as 



R + K -W^ 
-W qA + M 
-B D 



-B^ 






Uh 




fh 




o-h 









. 4>h . 








(12) 



where fh is the vector form of discretization of the linear form /(•). Since our goal is 
to obtain an efficient numerical scheme, we want to statically condense out the degree of 
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freedom associated with cr/j and cj)f^. Looking closely at the linear system (jl2p . we find that 
if the matrix D is diagonal, we can easily eliminate the degree of freedom corresponding 
to cTh and 0/^. This then leads to a formulation involving only one unknown Uh- Let 
{cpi,--- be the standard nodal finite element basis of Sh. We define a space 

spanned by the basis {fii,--- where the basis functions of and satisfy a 

condition of biorthogonality relation 



where n := dimM/i = dim Sh, ^ij is the Kronecker symbol, and Cj a positive scaling 
factor. This scaling factor Cj is chosen to be proportional to the area |supp(/5j[. In the 
following, we give these basis functions for linear simplicial finite elements in two and three 
dimensions. The d-parallelotope case is handled by using the tensor product construction 
of the one-dimensional basis functions presented in [35]. For the reference triangle T := 
{{x, y) : < X, < y, a; + y < 1}, we have 



where the basis functions /ii, jl2 and /is are associated with three vertices (0, 0), (1, 0) and 
(0, 1) of the reference triangle. For the reference tetrahedron T := {{x,y,z) : < a;,0 < 
y,0 < z,x + y + z < 1}, we have 

fli := A — 5x — 5y — 5z, (x2 := 5x — 1, and (x^ := by — 1, /i4 := bz — 1, 

where the basis functions /ii, /i2, As ^i^d ^4 associated with four vertices (0,0,0), (1,0,0), 
(0, 1,0) and (0,0, 1) of the reference tetrahedron. The global basis functions for the test 
space are constructed by glueing the local basis functions together following exactly the 
same procedure of constructing global finite element basis functions from the local ones. 
These global basis functions then satisfy the condition of biorthogonality (jl3p with global 
finite element basis functions. As these functions in are defined exactly in the same way 
as the finite element basis functions in Sh-, they satisfy supp^Uj = supp(/9j for i = 1, • • • , n. 
After statically condensing out variables (Th and (f)h-, we arrive at a reduced system 



Remark 1. Such biorthogonal basis functions are very popular in the context of mortar 
finite elements [23, 24, 35]. Construction of local basis functions of the space Mh satisfying 
all three AssumptionsU](i)^l^iii) as well as the biorthogonality condition (jlSh for different 
finite element spaces can be found in [26, 27, 35]. Working with nodal finite element basis 
functions based on Gauss-Lobatto quadrature nodes for rectangular or hexahedral triangu- 
lation, we have shown the construction of local basis functions of Mh satisfying all these 
assumptions for an arbitrary order finite element space [27]. Outside the context of mor- 
tar finite elements, these biorthogonal basis functions have first been used in the nearly 
incompressible elasticity in [25]. 




(13) 



fii := 3 — ix — 4y, /x2 := 4x — 1, and fl^ := 4y — 1 



((R + K) - (W'^D~^B + b'^D^^W) + B^D"^(aA + M)D"^B) Uh = fh- 
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3 An a priori error estimate 



Before proceeding to establish an a priori error estimate, we want to eliminate the gradient 
of the smoother and Lagrange multiplier 0^ from the saddle point problem (jlip . To 
this end, we introduce a quasi-projection operator: : L'^(Q) — Sh, which is defined as 

/ QhVfJ.hdx = / vnhdx, V G L^(r2), G Mf,. 
Jn Jn 

This type of operator is introduced in [31] to obtain the finite element interpolation of 
non-smooth functions satisfying boundary conditions, and is used in [7] in the context of 
mortar finite elements. The definition of Qh allows us to write the weak gradient as 

where operator is applied to vector Vu^ componentwise. We see that Qh is well- 
defined due to Assumptions [TJi) and mii) . Furthermore, the restriction of Qh to Sh is 
the identity. Hence Qh is a projection onto the space Sh- We note that Qh is not the 
orthogonal projection onto Sh but an oblique projection onto Sh- Oblique projectors are 
studied extensively in [20], and different expressions for the norm of oblique projections 
are provided in [32]. According to the biorthogonality relation between the basis functions 
of Sh and Mh, the action of operator Qh on a function v G L?'{^) can be written as 

QhV = Y.^ (14) 

which tells that the operator Qh is local in the sense to be given below, see also [1]. Let 
S{T') be the patch of an element T' & Th which is the interior of the closed set 

S{T') = [j{T e Th ■- dT ndT' ^ iD}. (15) 

Then Qh is local in the sense that for any v G L'^{Q), the value of QhV at any point 
in r G T/j only depends on the value of v in S{T) [1]. In the following, we will use a 
generic constant C, which will take different values at different places but will be always 
independent of the mesh-size h. The stability of Qh in L-^-norm is shown in the following 
lemma [23]. 

Lemma 1. Under Assumption{l^ii) , 

WQhvh^n) < C'lbllL2{n) for V G L^(0). (16) 
Proof- By Assumption [U^ii) 

fnlJ-hQhvdx fnUh^dx ^ ^ 

\\Qhv\\LHn)<P sup = /3 sup ^ < PWvh^in)- (17) 

□ 



7 



In the following, : L^{^) Sh will denote the L^-orthogonal projection onto Sh- It 
is well-known that the operator P/^ is stable in both L^- and i?^-norms. Using the stability 
of the operator Qh in the L'^-norm, and of the operator in the H^-novm, we can show 
that Qh is also stable in the H^-norm, see [24] for the locally quasi-uniform case. 

Lemma 2. Under Assumption{l^ii) , 

\Qhw\m{n) < C\w\min) for w e H'^{n). 
Proof. Using the L^-stability from Lemma [T] and the inverse inequality, we get for w € 

\Qhw\H^n) < \QhW - Phw\Hi(n) + \Phw\Hi(n) 

< C (^^\\Qh{w - Phw)\\L2(n) + \w\Hi{n)^ 

< C - Phw\\L2(n) + \Mm{n)^ < C\w\Hi{n)- 

□ 

The following lemma establishes the approximation property of operator for a 
function v G H'^{Q), see also [24]. 

Lemma 3. Under AssumptionU^ii) , there exists a constant C independent of the mesh- 
size h so that for v S H^~^^{^1), < s < 1, we have 



\\v - Qhv\\L2(n) < C/i^+''|u|j^s+i(n) 
\\v - Qhv\\m{n) < Ch^\v\Hs+i(n)- 

Proof. We start with a triangle inequality 

\\v - Qhv\\L2(n) < ll^' - Phvh^n) + \\PhV - Qhv\\L\n)- 
Since Qh acts as an identity on S^, we have 

\\v - Qhvh^in) < \\v - Phvh^n) + \\Qh{PhV - v)\\L2(^Qy 
Now we use the L^-stability of Qh from Lemma [1] to obtain 

\\v - Qhv\\L^(n) < C\\v - Phvh^^n)- 



(18) 



The first inequality of (jl8p follows by using the approximation property of the orthogonal 
projection Ph onto Sh, see [10]. The second inequality of (fTB|) is proved similarly using the 
stability of Qh in iif^-norm and the approximation property of the orthogonal projection 
Ph onto Sh- □ 
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Using the property of operator Qh, we can eliminate the degrees of freedom corre- 
sponding to so that our problem is to find 

Uh = min Jaivh), (19) 

where 

Ja{vh) = \\Pvhf + a\\\/{QhS/vh)\\l2^n) 

+ IIQ/iVuh — Vf/i||i2(f^) — 2 (Pv/i) z. 

In order to show that this problem has a unique solution, we define a P-inner product 
(•, •) p with 

{uh,Vh)p = {PuhfPvh + a Va-h : Vt/^ dx + / {<Th - Vuh) ■ [th - Vvh) dx, 

Jn Jn 

where ah = Qh{^Uh) and Th = Qh(yvh)- 

The following theorem shows that the P-inner product defines an inner product on the 
vector space Sh given by ([4]). 

Theorem 1. Let a > and Q C ft have at least three non-coUinear points for d = 2 and 
and four non-coplanar points for d = 3. Then the P-inner product defined above is an 
inner product on the vector space Sh ■ 

Proof. In order to show that the P-inner product is indeed an inner product, we have to 
prove the following properties of P-inner product: 

(1) {vh, Vh)p > 0, and {vh, Vh)p = if and only if Vh = 0, Vh ^ Sh, 

(2) {vh + Wh,Zh)p = {vh,Zh)p + {wh,Zh)p, Vh,Wh,Zh e Sh, 

(3) {vh, hz)p = b{vh, Zh)p, Vh e Sh, be R, 

(4) {vh,Wh)p = {wh,Vh)p, Vh,WheSh- 

It is trivial to show that the P-inner product satisfies the second, third and fourth proper- 
ties. It is also obvious that {vh,Vh)p > 0, and {vh,Vh)p = if Vh = 0. It remains to show 
that {vh,Vh)p = implies Vh = 0. We have {vh,Vh)p = WPvhW^ + a|| Vt/i|[|2(q) + \\Th - 
Vw/i||L2(n) with Th = Qhi^Vh)- Let {vh,Vh)p = 0. Then, WPvhW^ = 0, HVt/jHIj^^^ = 
and \\Th — ^Vh\\L2{n) = separately as they are all positive. Since Th is continu- 
ous, II VT/i||i2(-f^) = if and only if Th is a constant vector function in ^1. Similarly, 
||t/i — V?;/i||^2(Q) = implies that Vvh is also constant in fi, and thus Vh is a global linear 
function in Q. On the other hand, ||-P?^/i|| = implies that Vh is zero on Q C Cl, which 
contains at least three non-collinear points for d = 2 or four non-coplanar points for d = 3. 
Hence Vh is a global linear function which is zero at three non-collinear points for d = 2 
or four non-coplanar points for d = 3, and therefore, identically vanishes in 0. □ 
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The P-norm of an element Uh € Sh induced by the inner product {■,-)p is given by 
\\uh\\p = WPuhW^ + a||V(5/i(V'U/i)|[^2^^^ + WQhS/uh - Vii/ill^a^^^. Let the bihnear form 
a(-, ■) be defined as 

a{uh,Vh) = {PuhfPvh + a I Vah-VThdx+ / {an - Vuh) ■ {th - "^Vh) dy. 

Jci Jci 

with (Th = Qhi^Uh) and Th = QhC^Vh)- Since the biUnear form a(-, •) is symmetric, the 
minimization problem (|19p is equivalent to the variational problem of finding Uh G Sh such 
that [10, 16] 

a{uh,Vh) = f{vh), Vh G Sh- (20) 
Furthermore, the following corollary holds. 

Corollary 1. Under the assumptions of TheoremUl the variational problem ()20p admits 
a unique solution which depends continuously on the data. 

Proof. Since Uh,Vh £ Sh, it follows that \a{uh,Vh)\ < \\uh\\p\\vh\\p and \ f{vh)\ < C\\vh\\p. 
Moreover, using the definition of P-norm a{vh, Vh) = \\vh\\p, and thus a(-, •) is elliptic with 
repsect to the norm || ■ \\p. Hence our variational problem (j20p has a unique solution by 
Lax-Milgram Lemma [13, 16]. From the definition of the P-inner product, we have 

a{vh,Vh) = hhW'p, Vh G Sh, 
and thus, for the solution Uh G Sh, \\uh\\'p = f{uh)- D 



Remark 2. Using the unique solution Uh of the variational problem (120p . we have a unique 
solution {uh, (Th) 

The error estimate is obtained in the energy norm \\ ■ \\a induced by the bilinear form 
A{-, ■) defined as 



u,a)\\A := J||Pn||2 +a|cr|2 „ + _ Vn||2 {u,a) E C\^) x [H\n)Y. (21) 



iHi(n) ^ II" " "'iiL2(n) 

Theorem 2. Let u be the solution of continuous problem ([1]) with u S H^{^), a = Vu 
and (f) = a Act, and Uh be that of discrete problem ()20p with ah = Qh'^Uh. Then there 
exists a constant C > independent of the mesh-size h so that 

\\{u-Uh,o- - crh)\\A <C { inf \\{u - Wh, ct - eh)\\A + h\cf)\jjim) ] . 

\(«;ft,0,.)eKerBh / 

Proof. Here u, cr and cf) satisfy [14] 

A{iu,a),{v,T)) + B{(t>,{v,T)) = f{v), {v,t) € V, 
B{^I^,{u,<t)) =0, ^ G [L^m"- 

Let {wh, Oh) G KerBh so that {uh — Wh, (^h — ^h) G Ker Bh, and hence 

l|. n Ml ^ A{{Uh- Wh,(Th- Oh),{Vh,Th)) 

\\[Uh-Wh,(Th-dh)\\A< sup 7- ^ 

K,Th)GKcrBh \\{Vh,rh)\\A 
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Since A{{u -Uh,cT - ah), {vh,Th)) + B{cf), {vh,Th)) = for all {vh,Th) G KerBh, we have 

A{{uh - Wh, cTh - Oh), ivh,Th) = A{{u -Wh,cT- Oh), {vh,Th)) + A{{uh -u,ah- cr), {vh, 

= A{{u -Wh,a- - Oh), {vh,Th)) + -8(0, (vh, th))- 

Denoting the orthogonal projection of (p onto [M/^]'^ with respect to L'^-inner product by 
cj>f^, we have 

B{cf),{vh,Th)) = / {Th-Vvh)-i(f)-^h)dx<Ch\\Th-Vvh\\L^n)\((>\m(Q)- (22) 
Jn 

The result then follows by using the continuity of A{-, ■). □ 

Lemma 4. Let q G H^{Q) and q\s{T) ^■^ quadratic polynomial. Then 

QhVq = Vq on T. (23) 

Proof. If g is a quadratic polynomial in S{T), Vq is a vector with each component being 
a linear polynomial in S{T). Let h,--- Jng^j,-^ be indices of vertices in S{T) ordered in 
such a way that /i, • • • , In^ a-^e vertices of element T, and /nr+ij ' ' ' > ^ns^T) indices of 
rest vertices in S{T). Denoting the support of ipi by Si and using the expression of QhV 
from (I14p . we have 

""^ Jo fJ-i vdx 

Qhv\T = Y.^ V>k- (24) 



i=l 



Let u be a linear polynomial in S{T). Then v = Vi(pi. on S{T). We substitute this 

expression of v in (|24p and obtain 

QhV = y^vjifii, 
1=1 

which concludes that QhV = v in T. □ 
Lemma 5. Let Vh € Sh and u G H^(i}) with s > |. T/ien for all T Th 

WQh'^VhWL^^iT) < C\\VVh\\L^{S(T)), (25) 

and 

WQh'^huWL^^iT) < CW'^uWL^isiT)), (26) 
where Ih is the Lagrange interpolation operator. 

Proof. We note that it is sufficient to show 

\\Qhv\\Lo^(T) < C\\v\\loo(^s(t)), V e L°°(S(r)) 
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for (j25p . The estimate ()25p then follows by noting that Vvh S L°°{S{T)) and Q/i ^-cts on 
vector componentwise. The formula (|14p for QftW yields 



\Qhv\\ 



L°°(T) 



EL fii V dx 



l<i<n 
TCS, 



Ci 



L°°{T) 



where Si is the support of as before. Thus an application of the Cauchy-Schwarz 
inequality leads to 



/ jJLiV dx 




/ fiivdx 


< 


Jn 









< ll/^i||L2{S.)lhllL2(s^). 



So 



\\QhV\\L^{T) < 2^ ■ 

l<i<n 
TCS, 



V5i||L°°(T)- 



Since q is proportional to the area of Si, we estimate the L^-norm by the L°°-norm and 
use the quasi-uniformity assumption to obtain 

ll/^i||L2(5,;) 11^11^2(5.) < C Ci\\v\\L^(^Si)i 

where C is independent of the mesh-size. Moreover, using ||(/7j||j;^oo(7i) < 1, we have 

\\Qhv\\L^(T)<C ^ \\v\\loc^s,)- 

l<i<n 
TCS, 

Noting that element T is fixed and summation is restricted to those i's with T C Si, we 
have 

\\QhV\\L^{T) < C'll^llL°°(5(r))! 

where S(T) is as defined in (fTCj) . 

To obtain the estimate (|26p . we start with the mean value theorem as in [11] 



|V/fen||ioo(5.(y)) < 1 1 Vn 1 1^00(5.(7^)), 



and apply estimate ((25 



□ 



We note that Lemma [Hand [5] correspond to properties (Rl) and (R3) of the gradient 
recovery operator G stated in [1, Pages 72-73]. We show that operator has the same 
approximation property as operator G as stated in Theorem 4.1 of [1] by combining the 
arguments of [1, Theorem 4.1] and [11, Theorem 4.5]. 



Theorem 3. Let u G H^{S{T)) n H^{Q), T ^T^. Then we have 

\\Vu - Q/,V4n||i2(j.) < Ch^\\u\\Hz^s(T))^ 
where C > and is independent of h and u. 



(27) 
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Proof. As d e {2,3}, H^{S{T)) C C°{S{T)). We start with estimating the L^-norm by 
the L°°-norm and use Lemma [5] to obtain 

-II 

llVii - (5/iV4n||i2(r) < C/i2 ||Vn - Qh^hu\\L°°{T) 

< Chi (||Vn||L^(T) + ||Q/,V4n||ioo(r)) < C/i^ ||Vn||ioo(5.(T))- (28) 

Let (7 be a quadratic polynomial in S{T) to be selected later. Now applying Lemma H] in 
(j28]) . we obtain 

||Vn - Q/iV/hii||i2(7.) = ||V(n - - QhV4(n - 9)11^2(7^) 

< C/ii||V(n-g)||ioo(s(r)). (29) 

Applying Theorem 3.1.6 of [16], we can choose a quadratic polynomial q such that 

||V(u - ^)||l°°{s(t)) < C'/i"^/i^[u|H3(5(r))- (30) 

Finally, the result follows by combining ()29p with (I30p . □ 

Remark 3. Lemma ^is a polynomial reproduction property of operator Qf^. This re- 
production property is enough to prove (I30p when simplicial meshes or meshes of d- 
parallelotopes are used. However, in general quadrilateral and hexahedral meshes, the 
reproduction property of operator Qh is not enough to prove ([30j) . see [3]. 

Summing over all elements of 7/i, we obtain the following global estimate. 
Corollary 2. Let u e H^{n). Then we have 

||Vn - Q/iV//iu||j^2(Q) < Ch^\\u\\H3(n), (31) 
where C > and is independent of h. 

Theorem [2] shows that an optimal a priori estimate follows if there exist functions Vh 
and T/j with {vh,Th) G KerBh having the property 

||(n -Vh,cT - Th)\\A < Ch\\u\\Ha(n), 
which is guaranteed by the following theorem. 

Theorem 4. Under the assumptions of Theorem\^ there exists (f/i,T/j) € KerBh such 
that 

\\{u -Vh,(T - Th)\\A < Ch\\u\\H3(n)- (32) 

Proof. Let Vh be the Lagrange interpolation of u with respect to the mesh 7/j. Then it is 
well-known that 

\W-Vh\\m{n) <h^~''W\H^n), k = 0,l. (33) 

Moreover, 

\\P{u-Vh)\\' <h^\u\j,,^^y (34) 
13 



Let us recall the definition of the error in the energy norm 
\\{u -Vh,(T - Th)\\A = \J\\P{u - VhW + - '^?«l/fi(f7) + Ik -Th-Vu + yvh\\\2^p^^■ 

Let T/i = QhSvh so that {vh,Th) G KerBh- The approximation property of operator Qh 
given by Corollary [2] yields 

||Vn - (5hVf/i|li2(Q) < Ch^\\u\\H3{n)- (35) 
Hence, it suffices to show that 

Ik - '^hWn^n) < C'^ll'"llH3(n)- 
Since a = Vu and r/j = Qh^Vh, 

Ik ~ '''h\\m{n) < Ik — Qh(^\\m{n) + \\Qh(^ - Qh^vh\\m{n)- (36) 

The first term in the right-hand side of (j36|) has the correct approximation from Lemma 
[3l To estimate the second term, we use an inverse estimate 

C 

WQhO- - Qh'^VhWmin) < -j^WQhcr - Qh^Vh\\L2(^[^), (37) 

and apply the projection property and L^-stability of Qh to write 

C 

WQho- - Qh^vhWmin) < 'j^W'^^ ~ Qh^vhWi^in)- (38) 

Since Corollary [2] yields 

||Vn - QhVv/j||i2(Q) < Ch^\\u\\H3(n)^ (39) 

we have 

Ik - ''■/ill//i(n) < Ch\\u\\Hi{p.)- 

□ 

We combine the result of Theorems [2] and H] to get the final result. 

Theorem 5. Let u he the solution of continuous problem ^ with u G H^(yt), a = Vu 
and (f) = a Act, and Uh be that of discrete problem ()20p with = QhVuh. Then there 
exists a constant C > independent of the mesh-size h so that 

\\{u-Uh,(T -ah)\\A < Ch (||n||^3(f^) + \(t)\m{n)) ■ 

4 Conclusion 

We have presented a stabilized mixed finite element method for approximating thin plate 
splines in two and three dimensions. The mixed formulation introduces two additional 
vector variables - gradient of the smoother and Lagrange multiplier - as unknowns. In 
order to be able to eliminate these variables in an efficient way, we propose to use a pair of 
finite element bases satisfying a biorthogonality property for discretizing the gradient and 
the Lagrange multiplier. We have shown that the finite element approximation converges 
to the true solution of thin plate splines in an optimal way. 
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